This notebook draws from previous notebooks on relative cumulative frequency analysis and PCA analysis of metabolic cage data.
RMR data before addition of the wheels was analyzed and presented previously.
These datasets are from RMR data collected for mice after introduction to the metabolic cages (days 0-6, see image).
A total of 24 mice were run over the length of 2 experiments.
Experiment 1: The first set of mice, run in February
| Start of Met Cage run |
Metabolic Cage # |
Ear Tag | DOB | Gender | Genotype | Body Mass on BMR day (g) |
|---|---|---|---|---|---|---|
| 2/4/2020 | 1 | 1049 | 10/08/19 | M | WT | 27.4 |
| 2/4/2020 | 2 | 1050 | 10/08/19 | M | L-Bmal1-KO | 26.6 |
| 2/4/2020 | 3 | 1053 | 10/07/19 | M | WT | 31.5 |
| 2/4/2020 | 4 | 1054 | 10/07/19 | M | WT | 30.5 |
| 2/4/2020 | 5 | 1055 | 10/07/19 | M | L-Bmal1-KO | 27.5 |
| 2/4/2020 | 6 | 1109 | 10/11/19 | M | WT | 28.3 |
| 2/4/2020 | 9 | 1110 | 10/11/19 | M | WT | 27.9 |
| 2/4/2020 | 10 | 1111 | 10/11/19 | M | L-Bmal1-KO | 26.7 |
| 2/4/2020 | 11 | 1112 | 10/11/19 | M | WT | 28.5 |
| 2/4/2020 | 12 | 1214 | 10/21/19 | M | L-Bmal1-KO | 25.3 |
| 2/4/2020 | 13 | 1215 | 10/21/19 | M | WT | 27.9 |
| 2/4/2020 | 14 | 1216 | 10/21/19 | M | L-Bmal1-KO | 27.6 |
Experiment 2: The second set of mice, run in March
| Start of Met Cage run |
Metabolic Cage # |
Ear Tag | DOB | Gender | Genotype | Body Mass on BMR day (g) |
|---|---|---|---|---|---|---|
| 3/10/2020 | 1 | 1555 | 11/18/19 | M | WT | 26.4 |
| 3/10/2020 | 2 | 1692 | 12/06/19 | M | WT | 25.6 |
| 3/10/2020 | 3 | 1693 | 12/06/19 | M | WT | 24.7 |
| 3/10/2020 | 4 | 1695 | 12/06/19 | M | WT | 25.6 |
| 3/10/2020 | 5 | 1699 | 12/06/19 | M | WT | 28.7 |
| 3/10/2020 | 6 | 1556 | 11/18/19 | M | L-Bmal1-KO | 26.2 |
| 3/10/2020 | 9 | 1557 | 11/18/19 | M | L-Bmal1-KO | 27.1 |
| 3/10/2020 | 10 | 1558 | 11/18/19 | M | L-Bmal1-KO | 26.9 |
| 3/10/2020 | 11 | 1559 | 11/18/19 | M | L-Bmal1-KO | 24.7 |
| 3/10/2020 | 12 | 1691 | 12/06/19 | M | L-Bmal1-KO | 28 |
| 3/10/2020 | 13 | 1694 | 12/06/19 | M | L-Bmal1-KO | 26.4 |
| 3/10/2020 | 14 | 1700 | 12/06/19 | M | L-Bmal1-KO | 27.7 |
Data generated using a python script, called via shell command in R. A conda environment is required (see PCA notes ).
Bmal.py
#!/usr/bin/env python
"""
A python script that uses the mousetrapper library to prepare WT vs KO mice data for import to R.
Must be run from same directory.
"""
import numpy as np
import pandas as pd
import lib.mousetrapper as mt
__author__ = "Sumeed Yoyo Manzoor"
__email__ = "smanzoor@uchicago.edu"
__version__ = "0.0.1"
__experiment__ = "12 WT vs 12 L-Bmal1-KO mice"
datasets = "../datasets/"
data = "../data/"
outfolder = data + "R/"
experiment1_RMR = mt.Experiment.experiment(
info="""First run with revised protocol, WT vs L-Bmal1-KO""",
RT_data=datasets + "Bmal_WT_vs_KO_SPF_RMR_Wheel_2020-02-10_rt.csv",
macro13_data=datasets + "bmal_wt_vs_ko_spf_rmr_wheel_fixedwheel_2020-02-10_macro_One-Click Macro V2.32.2-slice3-VOC.mac_1.csv",
t_ctrl={"SPF WT":[1,3,4,6,9,11,13]},
t_exp={"SPF L-Bmal1-KO":[2,5,10,12,14]},
groups={
"SPF WT":[1,3,4,6,9,11,13],
"SPF L-Bmal1-KO":[2,5,10,12,14]
},
meta=["WT","L-Bmal1-KO","WT","WT","L-Bmal1-KO","WT","WT","L-Bmal1-KO","WT","L-Bmal1-KO","WT","L-Bmal1-KO"],
covar=[27.4,26.6,31.5,30.5,27.5,28.3,27.9,26.7,28.5,25.3,27.9,27.6]
)
experiment2_RMR = mt.Experiment.experiment(
info="""Second run with revised protocol, WT vs L-Bmal1-KO""",
RT_data=datasets + "Bmal_SPF_KF_RMR_Wheel_2020-03-16_16_31_rt.csv",
macro13_data=datasets + "bmal_spf_kf_rmr_wheel_2020-03-16_16_31_fixedwheel_macro_One-Click Macro V2.32.2-slice3-VOC.mac_1.csv",
t_ctrl={"SPF WT":np.arange(1,6).tolist()},
t_exp={"SPF L-Bmal1-KO":[6] + np.arange(9,15).tolist()},
groups={
"SPF WT":np.arange(1,6).tolist(),
"SPF L-Bmal1-KO":[6] + np.arange(9,15).tolist()
},
meta=np.array(["WT","L-Bmal1-KO"]).repeat([5,7]).tolist(),
covar=[26.4,25.6,24.7,25.6,28.7,26.2,27.1,26.9,24.7,28,26.4,27.7]
)
experiment1_RMR.macro13_lightcycles(days=2)
experiment2_RMR.macro13_lightcycles(days=2)
experiment1_RMR.macro13_trimdays(days=2, timeformat=False)
experiment2_RMR.macro13_trimdays(days=2, timeformat=False)
experiment1_RMR.macro13_dark.to_csv(outfolder + "1_dark.csv", index=False)
experiment1_RMR.macro13_light.to_csv(outfolder + "1_light.csv", index=False)
experiment1_RMR.macro13_trimmeddays.to_csv(outfolder + "1_total.csv", index=False)
pd.DataFrame({"covar":experiment1_RMR.covar}).to_csv(outfolder + "1_covar.csv", index=False)
experiment2_RMR.macro13_dark.to_csv(outfolder + "2_dark.csv", index=False)
experiment2_RMR.macro13_light.to_csv(outfolder + "2_light.csv", index=False)
experiment2_RMR.macro13_trimmeddays.to_csv(outfolder + "2_total.csv", index=False)
pd.DataFrame({"covar":experiment2_RMR.covar}).to_csv(outfolder + "2_covar.csv", index=False)
shell("conda activate met_cages_development && cd py && python Bmal_wheel.py")
reqpkg <- c("magrittr","dplyr","ggplot2","ggpubr","DT","reshape2")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "magrittr"
## [1] '1.5'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
## [1] "DT"
## [1] '0.13'
## [1] "reshape2"
## [1] '1.4.3'
theme_set(theme_pubr())
select <- dplyr::select
path <- "data/R/"
exp1 <- read.csv2(paste0(path, "1_total.csv"), sep = ",") %>%
set_colnames(paste0("exp1_",colnames(.)))
exp1_dark <- read.csv2(paste0(path, "1_dark.csv"), sep = ",") %>%
set_colnames(paste0("exp1_",colnames(.)))
exp1_light <- read.csv2(paste0(path, "1_light.csv"), sep = ",") %>%
set_colnames(paste0("exp1_",colnames(.)))
exp1_covar <- read.csv2(paste0(path, "1_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
exp2 <- read.csv2(paste0(path, "2_total.csv"), sep = ",") %>%
set_colnames(paste0("exp2_",colnames(.)))
exp2_dark <- read.csv2(paste0(path, "2_dark.csv"), sep = ",") %>%
set_colnames(paste0("exp2_",colnames(.)))
exp2_light <- read.csv2(paste0(path, "2_light.csv"), sep = ",") %>%
set_colnames(paste0("exp2_",colnames(.)))
exp2_covar <- read.csv2(paste0(path, "2_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
covar <- c(exp1_covar[c(1,3,4,6,7,9,11),], exp2_covar[1:5,], exp1_covar[c(2,5,8,10,12),], exp2_covar[6:12,])
total2 <- exp1 %>%
bind_cols(exp2) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
light <- exp1_light %>%
bind_cols(exp2_light) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
dark <- exp1_dark %>%
bind_cols(exp2_dark) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
total <- bind_rows(light, dark)
Select a tab above to view
light %>% select(contains("VO2")) %>% datatable(filter="bottom")
Note some data has errors (greater than 1 or below 0.5). These were likely measurement errors.
light %>% select(contains("RQ")) %>% datatable(filter="bottom")
light %>% select(contains("BodyMass")) %>% datatable(filter="bottom")
light %>% select(contains("Food")) %>% datatable(filter="bottom")
light %>% select(contains("Water")) %>% datatable(filter="bottom")
light %>% select(contains("PedMeters")) %>% datatable(filter="bottom")
light %>% select(contains("WheelMeters")) %>% datatable(filter="bottom")
Select a tab above to view
dark %>% select(contains("VO2")) %>% datatable(filter="bottom")
dark %>% select(contains("RQ")) %>% datatable(filter="bottom")
dark %>% select(contains("BodyMass")) %>% datatable(filter="bottom")
dark %>% select(contains("Food")) %>% datatable(filter="bottom")
dark %>% select(contains("Water")) %>% datatable(filter="bottom")
dark %>% select(contains("PedMeters")) %>% datatable(filter="bottom")
dark %>% select(contains("WheelMeters")) %>% datatable(filter="bottom")
select_WT <- function(df) {
d1 <- df %>% select(contains("exp1"))
for (name in colnames(d1)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(2,5,10,12,14)) {
d1 %<>% select(-all_of(name))
}
}
d2 <- df %>% select(contains("exp2"))
for (name in colnames(d2)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) > 5) {
d2 %<>% select(-all_of(name))
}
}
d <- bind_cols(d1, d2)
return(d)
}
select_KO <- function(df) {
d1 <- df %>% select(contains("exp1"))
for (name in colnames(d1)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,3,4,6,9,11,13)) {
d1 %<>% select(-all_of(name))
}
}
d2 <- df %>% select(contains("exp2"))
for (name in colnames(d2)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) <= 5) {
d2 %<>% select(-all_of(name))
}
}
d <- bind_cols(d1, d2)
return(d)
}
graph_groups <- function(data, channel, title="", type="m") {
df <- data %>% select(contains(channel))
df.observation <- 1:nrow(df)
c <- ncol(df)
WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt() %>% select(-variable)
KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt() %>% select(-variable)
df <- bind_cols(WT, KO) %>% set_colnames(c("WT","KO")) %>% melt()
df$observation <- rep(df.observation, c)
if (type=="m") {
df %>% group_by(variable, observation) %>%
summarise(value=mean(value)) %>%
ggplot(aes(x=observation, y=value, color=variable)) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_line(size=1) +
scale_color_manual(values = c("#103EFB", "#FC2A1C")) +
ggtitle(title) +
ylab(channel)
} else if (type=="s") {
ggplot(df, aes(x=observation, y=value, color=variable)) +
ggtitle(title) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
scale_color_manual(values = c("#103EFB", "#FC2A1C")) +
stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", alpha = 0.7) +
ylab(channel)
}
}
graph_mice <- function(data, channel, title="") {
df <- data %>% select(contains(channel))
df.observation <- 1:nrow(df)
c <- ncol(df)
WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt()
KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt()
df <- bind_rows(WT, KO) %>% set_colnames(c("mouse", "channel"))
df$observation <- rep(df.observation, c)
ggplot(df, aes(x=observation, y=channel, color=mouse)) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_line() +
ylab(channel)
}
There are 3 graphs per variable, a graph with the average of the mice in each group, a graph with the confidence limits for the population means (basically standard error), and a graph with each individual mouse plotted.
Select a tab above to view
graph_groups(total2, "VO2", title="VO2 by group", type = "m")
graph_groups(total2, "VO2", title="VO2 by group", type = "s")
graph_mice(total2, "VO2", title="VO2 by mouse")
total2 %>% filter_at(vars(contains("RQ")), all_vars(. < 1 & . > 0.5)) %>% graph_groups("RQ", title="RQ by group", type = "m")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
total2 %>% filter_at(vars(contains("RQ")), all_vars(. < 1 & . > 0.5)) %>% graph_groups("RQ", title="RQ by group", type = "s")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
total2 %>% filter_at(vars(contains("RQ")), all_vars(. < 1 & . > 0.5)) %>% graph_mice("RQ", title="RQ by mouse")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
graph_groups(total2, "BodyMass", title="BodyMass by group", type = "m")
graph_groups(total2, "BodyMass", title="BodyMass by group", type = "s")
graph_mice(total2, "BodyMass", title="BodyMass by mouse")
graph_groups(total2, "Food", title="Food by group", type = "m")
graph_groups(total2, "Food", title="Food by group", type = "s")
graph_mice(total2, "Food", title="Food by mouse")
graph_groups(total2, "Water", title="Water by group", type = "m")
graph_groups(total2, "Water", title="Water by group", type = "s")
graph_mice(total2, "Water", title="Water by mouse")
graph_groups(total2, "PedMeters", title="PedMeters by group", type = "m")
graph_groups(total2, "PedMeters", title="PedMeters by group", type = "s")
graph_mice(total2, "PedMeters", title="PedMeters by mouse")
graph_groups(total2, "PedSpeed", title="PedSpeed by group", type = "m")
graph_groups(total2, "PedSpeed", title="PedSpeed by group", type = "s")
graph_mice(total2, "PedSpeed", title="PedSpeed by mouse")
graph_groups(total2, "AllMeters_R", title="AllMeters_R by group", type = "m")
graph_groups(total2, "AllMeters_R", title="AllMeters_R by group", type = "s")
graph_mice(total2, "AllMeters_R", title="AllMeters_R by mouse")
graph_groups(total2, "AllMeters_M", title="AllMeters_M by group", type = "m")
graph_groups(total2, "AllMeters_M", title="AllMeters_M by group", type = "s")
graph_mice(total2, "AllMeters_M", title="AllMeters_M by mouse")
graph_groups(total2, "WheelMeters", title="WheelMeters by group", type = "m")
graph_groups(total2, "WheelMeters", title="WheelMeters by group", type = "s")
graph_mice(total2, "WheelMeters", title="WheelMeters by mouse")
graph_groups(total2, "WheelSpeed", title="WheelSpeed by group", type = "m")
graph_groups(total2, "WheelSpeed", title="WheelSpeed by group", type = "s")
graph_mice(total2, "WheelSpeed", title="WheelSpeed by mouse")
graph_groups(total2, "ZBreak", title="ZBreak by group", type = "m")
graph_groups(total2, "ZBreak", title="ZBreak by group", type = "s")
graph_mice(total2, "ZBreak", title="ZBreak by mouse")
graph_groups(total2, "H2O", title="H2O by group", type = "m")
graph_groups(total2, "H2O", title="H2O by group", type = "s")
graph_mice(total2, "H2O", title="H2O by mouse")
graph_groups(total2, "VOC", title="VOC by group", type = "m")
graph_groups(total2, "VOC", title="VOC by group", type = "s")
graph_mice(total2, "VOC", title="VOC by mouse")
graph_groups(total2, "H2_", title="H2 by group", type = "m")
graph_groups(total2, "H2_", title="H2 by group", type = "s")
graph_mice(total2, "H2_", title="H2 by mouse")
The idea behind this was to apply a method to study large metabolic datasets quantitatively. Based on the Gordon paper, we can use relative cumulative frequency to find frequency of data binned by VO2 or RQ reading levels, then calculate EC50. We can use EC50 values for every mouse in a simple t-test or ANCOVA to compare groups. In other words, this can be thought of as a kind of weighted average.
Relative cumulative frequency
This is a method where data is first ordered by increasing value, then, for every increase along a set increment, a frequency is recorded. Finally, frequency is cumulated. This is all accomplished by the function rcf. In order to run this, however, the range of data needs to be checked to determine suitable range and increments. find_range can be run to check the range on a vector, and will return a recommended range and increment value to submit to rcf.
EC50
By fitting frequencies to intervals, a sigmoidal regression is found. This is dependent on a normal distribution — that is, the mouse should have less frequencies towards the extremes, and the most expression in the middle. For RQ and VO2, this is true, the mouse will not spend as much time on the extremes of the values. For other data types, the nature and distribution of the data should be considered before using in this analysis.
t-test and ANCOVA
For the Welch t-test, I run an F-test first (var.test), which should be not significant.
ANCOVA is run using the EC50 values as dependent variable and bodyweight (see datasets) as covariate. The ANCOVA function returns a slew of information related to the test. A linear regression is first fit to the dependent variable (DV) and covariate. This may be important (see the Speakman paper ) since it can help elucidate the validity of ANCOVA related to other models, but unfortunately without lean mass data the conclusions to be drawn are limited. Evidently, rarely does energy expenditure regress linearly with total mass. Next, ANOVA is run between DV and covariate. Hopefully, the covariate is not a significant predictor of DV. Then, a boxplot for DV vs independent variable (IV), in this case WT vs KO, and covariate vs IV are shown to check for extreme outliers and IV relationship. Then, basic stats on all variables are displayed, followed by Levene’s test, which should not have a significant result. Sanity checks are run using ANOVA of DV, ANOVA of covariate, then ANOVA with DV normalized to covariate (which should not be used). ANCOVA is finally run. The check_ANCOVA function will pull any significant (<0.05) p values from ANCOVA.
In the future, more complicated models, such as ANCOVA/linear models with more covariates, might be worth implementing.
reqpkg <- c("docstring","tibble","scales","purrr","nplr", "car","compute.es","effects","multcomp","pastecs","WRS2")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "docstring"
## [1] '1.0.0'
## [1] "tibble"
## [1] '2.1.3'
## [1] "scales"
## [1] '1.1.0'
## [1] "purrr"
## [1] '0.3.3'
## [1] "nplr"
## [1] '0.1.7'
## [1] "car"
## [1] '3.0.7'
## [1] "compute.es"
## [1] '0.2.4'
## [1] "effects"
## [1] '4.1.4'
## [1] "multcomp"
## [1] '1.4.12'
## [1] "pastecs"
## [1] '1.3.21'
## [1] "WRS2"
## [1] '1.0.0'
find_range <- function(data, channel, segments=20) {
data <- data %>% select(contains(channel)) %>% melt()
r <- range(data$value)
cut <- max(data$value)-min(data$value)
out <- list("range"=r, "cut"=cut/segments)
return(out)
}
rcf <- function(data, channel, id, min, max, b=0.1) {
#' Relative cumulative frequency function
#'
#' Takes a vector and returns a dataframe with relative cumulative frequency per break unit
#'
#' @param data numeric vector
#' @param channel character. Channel to label the rcf output.
#' @param b double. Number to cut vector.
#' Before running, check for valid break values to cut the data. One way is to divide range of data by 10 `(max(data)-min(data))/10`.
breaks <- seq(min, max, by=b)
duration.cut <- cut(data, breaks, right = F)
duration.freq <- table(duration.cut)
duration.cumfreq = cumsum(duration.freq) %>% prepend(0)
duration.cumrelfreq = duration.cumfreq/length(data)
out <- duration.cumrelfreq %>% as.data.frame()
out <- cbind(breaks, out) %>% set_colnames(c("breaks", channel)) %>% set_rownames(NULL)
out$id <- rep(id, nrow(out))
return(out)
}
EC50 <- function(x) {
out <- lapply(x, function(y) {
nplr(y$breaks, y[, 2], silent = T, useLog = FALSE)
}) %>%
sapply(function(z) {
getEstimates(z, 0.5)$x
})
out %<>% set_rownames(NULL)
return(out)
}
lm_eqn <- function(df, y, x){
m <- lm(y ~ x, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(unname(coef(m)[1]), digits = 2),
b = format(unname(coef(m)[2]), digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
lm_eqn_str <- function(df, y, x){
m <- lm(y ~ x, df);
eq <- sprintf("y = %s + %sx, r^2 = %s", format(unname(coef(m)[1]), digits = 2),format(unname(coef(m)[2]), digits = 2),format(summary(m)$r.squared, digits = 3))
return(eq)
}
remove_outliers <- function(x, na.rm = TRUE, ...) {
qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
H <- 1.5 * IQR(x, na.rm = na.rm)
y <- x
y[x < (qnt[1] - H)] <- NA
y[x > (qnt[2] + H)] <- NA
y
}
ANCOVA <- function(df) {
# Takes dataframe with column 1 as independent variable (e.g., WT), column 2 as dependent variable, and column 3 as covariate
cat("linear fit of dependent variable with covariates\n")
p1 <- ggplot(df, aes_string(x=names(df)[2], y=names(df)[3])) +
geom_smooth(method = "lm", se = FALSE) +
geom_point() +
# geom_text(x = mean(na.omit(df[[2]])), y = mean(na.omit(df[[3]])), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
geom_text(x = mean(na.omit(df[[2]])), y = (mean(na.omit(df[[3]])) + IQR(df[[3]], na.rm = T)), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
ggtitle("Linear model - DV and Covar") +
xlab("DV") + ylab("Covar")
print(p1)
m <- lm(df[[3]] ~ df[[2]], df)
summary(m) %>% print()
lm_eqn_str(df, df[[3]], df[[2]]) %>% cat()
plot1 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[2])) +
geom_boxplot() +
# ggtitle("DV vs IV") +
xlab("IV") + ylab("DV")
plot2 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[3])) +
geom_boxplot() +
# ggtitle("Covar vs IV") +
xlab("IV") + ylab("Covar")
p2 <- ggarrange(plot1, plot2,
labels = c("DV vs IV", "Covar vs IV"),
ncol = 2)
# grid.arrange(plot1, plot2, ncol=2)
print(p2)
cat("\n_____________________\n")
cat("Some stats on the variables\n")
cat("DV - IV:\n")
by(df[[2]], df[[1]], stat.desc) %>% print()
cat("Covar - IV:\n")
by(df[[3]], df[[1]], stat.desc) %>% print()
cat("levene's test\n")
leveneTest(df[[2]], df[[1]], center = median) %>% print()
cat("_____________________\n")
cat("ANOVA of DV with IV\n")
model <- aov(df[[2]] ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANOVA of covariate with IV, independent of DV\n")
model <- aov(df[[3]] ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANOVA of DV with IV, normalized to covariate\n")
c <- df[[2]]/df[[3]]
model <- aov(c ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANCOVA\n")
model <- aov(df[[2]] ~ df[[1]] + df[[3]], data = df)
summary(model) %>% print()
out <- Anova(model, type= "III")
print(out)
return(out)
}
check_ANCOVA <- function(res) {
if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] < 0.05) {
print(paste("Varies significantly with DV and Covar,"), res$`Pr(>F)`[2], res$`Pr(>F)`[3], sep=" ")
} else if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] > 0.05) {
print(paste("Varies significantly with DV,", "p =", res$`Pr(>F)`[2], sep =" "))
} else if (res$`Pr(>F)`[2] > 0.05 && res$`Pr(>F)`[3] < 0.05) {
print(paste("Varies significantly with Covar,", "p =", res$`Pr(>F)`[3], sep=" "))
}
}
Select a tab above to view
total.VO2 <- total %>% select(contains("VO2"))
total.VO2.WT <- total.VO2 %>% select_WT()
total.VO2.KO <- total.VO2 %>% select_KO()
rcf.WT <- lapply(total.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.001))
rcf.KO <- lapply(total.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.001))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})
# overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in complete cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 1.71539878847918"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 1.92269230171977"
Run F & t test
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.5506, num df = 11, denom df = 11, p-value = 0.4787
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.4463848 5.3863404
## sample estimates:
## ratio of variances
## 1.550607
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -1.7594, df = 21.02, p-value = 0.09306
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3565851 0.0297364
## sample estimates:
## mean of x mean of y
## 1.776285 1.939710
Run ANCOVA
res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()
## linear fit of dependent variable with covariates
## `geom_smooth()` using formula 'y ~ x'
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6789 -0.9095 -0.1433 0.7997 4.0970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.5857 2.7174 10.52 4.76e-10 ***
## df[[2]] -0.7256 1.4512 -0.50 0.622
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.654 on 22 degrees of freedom
## Multiple R-squared: 0.01124, Adjusted R-squared: -0.03371
## F-statistic: 0.25 on 1 and 22 DF, p-value: 0.622
##
## y = 29 + -0.73x, r^2 = 0.0112
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.58030227 2.30114728 0.72084501
## sum median mean SE.mean CI.mean.0.95 var
## 23.27651438 1.91601402 1.93970953 0.05816211 0.12801395 0.04059398
## std.dev coef.var
## 0.20147947 0.10387095
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.36140182 2.11187515 0.75047334
## sum median mean SE.mean CI.mean.0.95 var
## 21.31542198 1.67948358 1.77628517 0.07242541 0.15940726 0.06294529
## std.dev coef.var
## 0.25088899 0.14124365
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.792 0.3831
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.1602 0.16025 3.095 0.0924 .
## Residuals 22 1.1389 0.05177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0004118 0.0004118 5.162 0.0332 *
## Residuals 22 0.0017550 0.0000798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.1602 0.16025 2.955 0.100
## df[[3]] 1 0.0001 0.00007 0.001 0.971
## Residuals 21 1.1389 0.05423
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.27674 1 5.1030 0.03464 *
## df[[1]] 0.14572 1 2.6869 0.11607
## df[[3]] 0.00007 1 0.0013 0.97144
## Residuals 1.13886 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
light.VO2 <- light %>% select(contains("VO2"))
light.VO2.WT <- light.VO2 %>% select_WT()
light.VO2.KO <- light.VO2 %>% select_KO()
rcf.WT <- lapply(light.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.001))
rcf.KO <- lapply(light.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.001))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in light cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 1.38854018713538"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 1.44765116831737"
Run F & t test
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 0.44565, num df = 11, denom df = 11, p-value = 0.1959
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1282914 1.5480390
## sample estimates:
## ratio of variances
## 0.4456457
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -1.445, df = 19.18, p-value = 0.1646
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.24328399 0.04448241
## sample estimates:
## mean of x mean of y
## 1.373416 1.472816
Run ANCOVA
res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()
## linear fit of dependent variable with covariates
## `geom_smooth()` using formula 'y ~ x'
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4161 -1.0136 -0.0008 0.6401 4.0365
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.484 2.696 8.340 2.94e-08 ***
## df[[2]] 3.341 1.881 1.776 0.0896 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.556 on 22 degrees of freedom
## Multiple R-squared: 0.1254, Adjusted R-squared: 0.0856
## F-statistic: 3.153 on 1 and 22 DF, p-value: 0.08962
##
## y = 22 + 3.3x, r^2 = 0.125
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.21686523 1.92189357 0.70502834
## sum median mean SE.mean CI.mean.0.95 var
## 17.67379688 1.43491951 1.47281641 0.05721110 0.12592078 0.03927732
## std.dev coef.var
## 0.19818505 0.13456195
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.17590660 1.51816084 0.34225424
## sum median mean SE.mean CI.mean.0.95 var
## 16.48098736 1.39123252 1.37341561 0.03819224 0.08406055 0.01750377
## std.dev coef.var
## 0.13230180 0.09633049
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2399 0.6291
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0593 0.05928 2.088 0.163
## Residuals 22 0.6246 0.02839
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0001871 1.871e-04 6.57 0.0177 *
## Residuals 22 0.0006264 2.847e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0593 0.05928 2.644 0.1189
## df[[3]] 1 0.1537 0.15366 6.852 0.0161 *
## Residuals 21 0.4709 0.02243
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.00023 1 0.0101 0.92079
## df[[1]] 0.12721 1 5.6728 0.02677 *
## df[[3]] 0.15366 1 6.8522 0.01608 *
## Residuals 0.47093 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV and Covar,"
dark.VO2 <- dark %>% select(contains("VO2"))
dark.VO2.WT <- dark.VO2 %>% select_WT()
dark.VO2.KO <- dark.VO2 %>% select_KO()
rcf.WT <- lapply(dark.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.001))
rcf.KO <- lapply(dark.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.001))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in dark cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 2.24612143206034"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 2.59513184736794"
Run F & t test
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 0.93523, num df = 11, denom df = 11, p-value = 0.9136
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.2692306 3.2486944
## sample estimates:
## ratio of variances
## 0.9352262
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -2.1066, df = 21.975, p-value = 0.0468
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.80197969 -0.00625586
## sample estimates:
## mean of x mean of y
## 2.263029 2.667146
Run ANCOVA
res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()
## linear fit of dependent variable with covariates
## `geom_smooth()` using formula 'y ~ x'
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6998 -0.9217 -0.2470 0.9448 3.7247
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.1574 1.6797 17.358 2.52e-14 ***
## df[[2]] -0.7788 0.6682 -1.166 0.256
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.614 on 22 degrees of freedom
## Multiple R-squared: 0.05816, Adjusted R-squared: 0.01535
## F-statistic: 1.359 on 1 and 22 DF, p-value: 0.2563
##
## y = 29 + -0.78x, r^2 = 0.0582
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.0000000 0.0000000 0.0000000 1.8745230 3.8776946 2.0031716
## sum median mean SE.mean CI.mean.0.95 var
## 32.0057573 2.5865853 2.6671464 0.1378974 0.3035102 0.2281884
## std.dev coef.var
## 0.4776907 0.1791018
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.0000000 0.0000000 0.0000000 1.5858365 2.8022937 1.2164572
## sum median mean SE.mean CI.mean.0.95 var
## 27.1563440 2.3163711 2.2630287 0.1333566 0.2935159 0.2134077
## std.dev coef.var
## 0.4619608 0.2041338
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.5967 0.4481
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.980 0.9799 4.438 0.0468 *
## Residuals 22 4.858 0.2208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.001834 0.0018341 5.682 0.0262 *
## Residuals 22 0.007101 0.0003228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.980 0.9799 4.305 0.0505 .
## df[[3]] 1 0.078 0.0778 0.342 0.5649
## Residuals 21 4.780 0.2276
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 1.0262 1 4.5087 0.04577 *
## df[[1]] 0.7182 1 3.1554 0.09017 .
## df[[3]] 0.0778 1 0.3420 0.56492
## Residuals 4.7797 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
Note that some datapoints were removed due to measurement error. Data was filtered for values < 1 and > 0.5.
Select a tab above to view
total.RQ <- total %>% select(contains("RQ")) %>% filter_all(all_vars(. < 1 & . > 0.5))
total.RQ.WT <- total.RQ %>% select_WT()
total.RQ.KO <- total.RQ %>% select_KO()
rcf.WT <- lapply(total.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.0005))
rcf.KO <- lapply(total.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.0005))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in complete cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 0.815989968208129"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 0.820796165978486"
Run F & t test
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.7748, num df = 11, denom df = 11, p-value = 0.3555
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.510926 6.165133
## sample estimates:
## ratio of variances
## 1.774803
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -0.72408, df = 20.409, p-value = 0.4772
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.02871452 0.01390229
## sample estimates:
## mean of x mean of y
## 0.8160276 0.8234337
Run ANCOVA
res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()
## linear fit of dependent variable with covariates
## `geom_smooth()` using formula 'y ~ x'
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6588 -0.9693 0.2279 0.8568 4.1602
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.90 11.30 1.584 0.127
## df[[2]] 11.39 13.78 0.827 0.417
##
## Residual standard error: 1.638 on 22 degrees of freedom
## Multiple R-squared: 0.03015, Adjusted R-squared: -0.01394
## F-statistic: 0.6838 on 1 and 22 DF, p-value: 0.4172
##
## y = 18 + 11x, r^2 = 0.0301
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 1.200000e+01 0.000000e+00 0.000000e+00 7.946162e-01 8.641973e-01 6.958113e-02
## sum median mean SE.mean CI.mean.0.95 var
## 9.881205e+00 8.237821e-01 8.234337e-01 6.140242e-03 1.351458e-02 4.524308e-04
## std.dev coef.var
## 2.127042e-02 2.583137e-02
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 1.200000e+01 0.000000e+00 0.000000e+00 7.683244e-01 8.531133e-01 8.478887e-02
## sum median mean SE.mean CI.mean.0.95 var
## 9.792331e+00 8.192496e-01 8.160276e-01 8.180137e-03 1.800436e-02 8.029757e-04
## std.dev coef.var
## 2.833683e-02 3.472533e-02
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.7314 0.4016
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000329 0.0003291 0.524 0.477
## Residuals 22 0.013809 0.0006277
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 1.087e-05 1.087e-05 3.464 0.0761 .
## Residuals 22 6.902e-05 3.137e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000329 0.0003291 0.531 0.474
## df[[3]] 1 0.000782 0.0007823 1.261 0.274
## Residuals 21 0.013027 0.0006203
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.039605 1 63.8444 8.398e-08 ***
## df[[1]] 0.000685 1 1.1046 0.3052
## df[[3]] 0.000782 1 1.2611 0.2741
## Residuals 0.013027 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
light.RQ <- light %>% select(contains("RQ")) %>% filter_all(all_vars(. < 1 & . > 0.5))
light.RQ.WT <- light.RQ %>% select_WT()
light.RQ.KO <- light.RQ %>% select_KO()
rcf.WT <- lapply(light.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.0005))
rcf.KO <- lapply(light.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.0005))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in light cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 0.772440841371076"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 0.780156358993069"
Run F & t test
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 3.2866, num df = 11, denom df = 11, p-value = 0.06047
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.9461303 11.4165622
## sample estimates:
## ratio of variances
## 3.286572
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -0.32483, df = 17.127, p-value = 0.7492
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.03087977 0.02263589
## sample estimates:
## mean of x mean of y
## 0.7740705 0.7781924
Run ANCOVA
res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()
## linear fit of dependent variable with covariates
## `geom_smooth()` using formula 'y ~ x'
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2970 -1.1931 0.2283 1.0840 3.7131
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.378 8.168 1.393 0.1775
## df[[2]] 20.434 10.516 1.943 0.0649 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.537 on 22 degrees of freedom
## Multiple R-squared: 0.1465, Adjusted R-squared: 0.1077
## F-statistic: 3.776 on 1 and 22 DF, p-value: 0.0649
##
## y = 11 + 20x, r^2 = 0.146
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.000000000 0.000000000 0.000000000 0.734585496 0.805507815 0.070922319
## sum median mean SE.mean CI.mean.0.95 var
## 9.338309179 0.777979674 0.778192432 0.006129090 0.013490037 0.000450789
## std.dev coef.var
## 0.021231791 0.027283472
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 0.71307952 0.85084419 0.13776467
## sum median mean SE.mean CI.mean.0.95 var
## 9.28884594 0.77121792 0.77407049 0.01111137 0.02445596 0.00148155
## std.dev coef.var
## 0.03849091 0.04972533
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.4256 0.1336
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000102 0.0001019 0.106 0.748
## Residuals 22 0.021256 0.0009662
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 8.420e-06 8.418e-06 3.459 0.0763 .
## Residuals 22 5.354e-05 2.434e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000102 0.000102 0.123 0.7289
## df[[3]] 1 0.003907 0.003907 4.730 0.0412 *
## Residuals 21 0.017348 0.000826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.0231390 1 28.0095 3.015e-05 ***
## df[[1]] 0.0008806 1 1.0659 0.31361
## df[[3]] 0.0039074 1 4.7299 0.04122 *
## Residuals 0.0173483 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with Covar, p = 0.0412152943838807"
dark.RQ <- dark %>% select(contains("RQ")) %>% filter_all(all_vars(. < 1 & . > 0.5))
dark.RQ.WT <- dark.RQ %>% select_WT()
dark.RQ.KO <- dark.RQ %>% select_KO()
rcf.WT <- lapply(dark.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.0005))
rcf.KO <- lapply(dark.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.0005))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in dark cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 0.844890762938068"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 0.850790371295462"
Run F & t test
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.1212, num df = 11, denom df = 11, p-value = 0.8529
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3227626 3.8946419
## sample estimates:
## ratio of variances
## 1.12118
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -0.8451, df = 21.928, p-value = 0.4072
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.03766279 0.01585748
## sample estimates:
## mean of x mean of y
## 0.8454874 0.8563901
Run ANCOVA
res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()
## linear fit of dependent variable with covariates
## `geom_smooth()` using formula 'y ~ x'
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5600 -0.8905 0.0512 0.7171 4.2693
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.302 9.403 2.797 0.0105 *
## df[[2]] 1.100 11.042 0.100 0.9216
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.663 on 22 degrees of freedom
## Multiple R-squared: 0.0004505, Adjusted R-squared: -0.04498
## F-statistic: 0.009915 on 1 and 22 DF, p-value: 0.9216
##
## y = 26 + 1.1x, r^2 = 0.00045
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.000000000 0.000000000 0.000000000 0.828795893 0.937811179 0.109015286
## sum median mean SE.mean CI.mean.0.95 var
## 10.276680821 0.847528600 0.856390068 0.008857991 0.019496306 0.000941568
## std.dev coef.var
## 0.030684981 0.035830612
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.000000000 0.000000000 0.000000000 0.801022574 0.920923635 0.119901060
## sum median mean SE.mean CI.mean.0.95 var
## 10.145848961 0.847581860 0.845487413 0.009379353 0.020643818 0.001055667
## std.dev coef.var
## 0.032491033 0.038428760
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.158 0.6949
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000713 0.0007132 0.714 0.407
## Residuals 22 0.021970 0.0009986
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 1.357e-05 1.357e-05 3.234 0.0859 .
## Residuals 22 9.231e-05 4.196e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000713 0.0007132 0.687 0.417
## df[[3]] 1 0.000155 0.0001551 0.149 0.703
## Residuals 21 0.021815 0.0010388
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.049979 1 48.1129 7.487e-07 ***
## df[[1]] 0.000858 1 0.8260 0.3737
## df[[3]] 0.000155 1 0.1493 0.7031
## Residuals 0.021815 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
reqpkg <- c("ggfortify","rgl","lfda")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "ggfortify"
## [1] '0.4.10'
## [1] "rgl"
## [1] '0.100.54'
## [1] "lfda"
## [1] '1.1.3'
knitr::knit_hooks$set(webgl = hook_webgl)
pca_mutate <- function(data, channel) {
df <- data %>% select(contains(channel))
WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt() %>% select(-variable)
KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt() %>% select(-variable)
df <- bind_cols(WT, KO) %>% set_colnames(c("WT","KO")) %>% melt()
return(df)
}
plot_pca <- function(pca, data, title="") {
autoplot(pca, data = na.omit(data), colour = "id", loadings = TRUE, loadings.colour = "black", loadings.label = TRUE, loadings.label.colour = "black", loadings.label.size = 5) + ggtitle(title)
}
Small filtering step: during the day, there were points where the gas analyzers failed, leaving VO2 and RQ values unbelievably low. The rows contianing data for those points were simply discarded.
light <- light %>% slice(-320:-325)
VO2_light <- pca_mutate(light, "VO2")
RQ_light <- light %>% pca_mutate("RQ")
Food_light <- pca_mutate(light, "Food")
Water_light <- pca_mutate(light, "Water")
BM_light <- pca_mutate(light, "BodyMass")
PedMeters_light <- pca_mutate(light, "PedMeters_R")
WheelMeters_light <- pca_mutate(light, "WheelMeters")
df_light <- RQ_light %>% set_colnames(c("id","RQ")) %>% add_column(VO2=VO2_light$value) %>% add_column(FoodIn=Food_light$value) %>% add_column(WaterIn=Water_light$value) %>% add_column(BodyMass=BM_light$value) %>% add_column(PedMeters=PedMeters_light$value) %>% add_column(WheelMeters=WheelMeters_light$value)
PCA of most variables
df_light %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title = "PCA of most variables, light cycle")
PCA without BodyMass
df_light %>% select(-id,-BodyMass) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA of most variables besides BM, light cycle")
VO2 vs RQ
df_light %>% select(VO2,RQ) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA VO2 and RQ, light cycle")
Food vs Water
df_light %>% select(FoodIn,WaterIn) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA Food and Water, light cycle")
Including extra variables (AllMeters, ZBreaks)
AllMeters_light <- pca_mutate(light, "AllMeters_R")
ZB_light <- pca_mutate(light, "ZBreak")
df_light %<>% add_column(AllMeters=AllMeters_light$value) %>% add_column(ZBreak=ZB_light$value)
df_light %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA of all variables, light cycle")
VO2_dark <- pca_mutate(dark, "VO2")
RQ_dark <- dark %>% pca_mutate("RQ")
Food_dark <- pca_mutate(dark, "Food")
Water_dark <- pca_mutate(dark, "Water")
BM_dark <- pca_mutate(dark, "BodyMass")
PedMeters_dark <- pca_mutate(dark, "PedMeters_R")
WheelMeters_dark <- pca_mutate(dark, "WheelMeters")
df_dark <- RQ_dark %>% set_colnames(c("id","RQ")) %>% add_column(VO2=VO2_dark$value) %>% add_column(FoodIn=Food_dark$value) %>% add_column(WaterIn=Water_dark$value) %>% add_column(BodyMass=BM_dark$value) %>% add_column(PedMeters=PedMeters_dark$value) %>% add_column(WheelMeters=WheelMeters_dark$value)
PCA of most variables
df_dark %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, "PCA of most variables, dark cycle")
PCA without BodyMass
df_dark %>% select(-id,-BodyMass) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA of most variables besides BM, dark cycle")
VO2 vs RQ
df_dark %>% select(VO2,RQ) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA VO2 and RQ, dark cycle")
Food vs Water
df_dark %>% select(FoodIn,WaterIn) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA Food and Water, dark cycle")
Including extra variables (AllMeters, ZBreaks)
AllMeters_dark <- pca_mutate(dark, "AllMeters_R")
ZB_dark <- pca_mutate(dark, "ZBreak")
df_dark %<>% add_column(AllMeters=AllMeters_dark$value) %>% add_column(ZBreak=ZB_dark$value)
df_dark %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA of all variables, dark cycle")
I attempted to model and plot ldfa. LDFA should give similar results to PCA, with the advantage of being sensitive to multimodality of data, i.e. multiple effects driving an observation, as well as sensitive to distinctiveness in the data that might be lost in reducing dimensions in PCA. In other words, a learning model like LDFA might account for more relationships than a principal component.
More information:
model <- lfda(df_light[-1], df_light[, 1], r=3, metric = "plain")
autoplot(model, data=df_light, frame = T, colour = "id")
plot(model, labels = df_light[, 1])
You must enable Javascript to view this page properly.
model <- lfda(df_dark[-1], df_dark[, 1], r=3, metric = "plain")
autoplot(model, data=df_dark, frame = T, colour = "id")
plot(model, labels = df_dark[, 1])
You must enable Javascript to view this page properly.
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lfda_1.1.3 rgl_0.100.54 ggfortify_0.4.10 WRS2_1.0-0
## [5] pastecs_1.3.21 multcomp_1.4-12 TH.data_1.0-10 MASS_7.3-51.5
## [9] survival_3.1-11 mvtnorm_1.1-0 effects_4.1-4 compute.es_0.2-4
## [13] car_3.0-7 carData_3.0-3 nplr_0.1-7 purrr_0.3.3
## [17] scales_1.1.0 tibble_2.1.3 docstring_1.0.0 reshape2_1.4.3
## [21] DT_0.13 ggpubr_0.2.5 ggplot2_3.3.0 dplyr_0.8.5
## [25] magrittr_1.5
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_1.4-1 ggsignif_0.6.0
## [4] rio_0.5.16 htmlTable_1.13.3 markdown_1.1
## [7] base64enc_0.1-3 mc2d_0.1-18 rstudioapi_0.11
## [10] roxygen2_7.1.0 farver_2.0.3 RSpectra_0.16-0
## [13] xml2_1.2.5 codetools_0.2-16 splines_3.6.3
## [16] knitr_1.28 Formula_1.2-3 jsonlite_1.6.1
## [19] nloptr_1.2.2.1 cluster_2.1.0 png_0.1-7
## [22] shiny_1.4.0.2 compiler_3.6.3 backports_1.1.5
## [25] fastmap_1.0.1 assertthat_0.2.1 Matrix_1.2-18
## [28] survey_3.37 later_1.0.0 acepack_1.4.1
## [31] htmltools_0.4.0 tools_3.6.3 gtable_0.3.0
## [34] glue_1.3.2 Rcpp_1.0.4.6 cellranger_1.1.0
## [37] vctrs_0.2.4 nlme_3.1-144 crosstalk_1.1.0.1
## [40] xfun_0.12 stringr_1.4.0 openxlsx_4.1.4
## [43] lme4_1.1-21 miniUI_0.1.1.1 mime_0.9
## [46] lifecycle_0.2.0 klippy_0.0.0.9500 zoo_1.8-7
## [49] hms_0.5.3 promises_1.1.0 sandwich_2.5-1
## [52] RColorBrewer_1.1-2 yaml_2.2.1 curl_4.3
## [55] gridExtra_2.3 rpart_4.1-15 reshape_0.8.8
## [58] latticeExtra_0.6-29 stringi_1.4.6 checkmate_2.0.0
## [61] manipulateWidget_0.10.1 boot_1.3-24 zip_2.0.4
## [64] rlang_0.4.5 pkgconfig_2.0.3 evaluate_0.14
## [67] lattice_0.20-40 htmlwidgets_1.5.1 labeling_0.3
## [70] cowplot_1.0.0 tidyselect_1.0.0 plyr_1.8.6
## [73] R6_2.4.1 Hmisc_4.4-0 DBI_1.1.0
## [76] pillar_1.4.3 haven_2.2.0 foreign_0.8-75
## [79] withr_2.1.2 mgcv_1.8-31 abind_1.4-5
## [82] nnet_7.3-13 crayon_1.3.4 rARPACK_0.11-0
## [85] rmarkdown_2.1 jpeg_0.1-8.1 grid_3.6.3
## [88] readxl_1.3.1 data.table_1.12.8 forcats_0.5.0
## [91] webshot_0.5.2 digest_0.6.25 xtable_1.8-4
## [94] tidyr_1.0.2 httpuv_1.5.2 munsell_0.5.0
## [97] mitools_2.4